# include "scene.h"
# include "discr.h"
# include <stdio.h>
# include <errno.h>
# include <string.h>
# include <math.h>
// # include <ctype.h>

extern double kvisc(double T);
extern double dens(double T);
extern double thm_cond(double T);

extern short print_result(unsigned int t_cnt, const struct grid *g)
  {
   int _errno;
   unsigned short i, j, k, nr, nth, nz, isec;
   char resfile[MAX_FNAME_LEN+1];
   FILE *ptr;
   double ***u, ***v, ***w, ***p, *r, *th, *z, _u, _v, _w, v_total, vort_th, div;
   double T_wall, stress_r, stress_th;
#  ifdef SOLVE_ENERGY
   double ***T, h, A_total, T_wall_av, h_av, q_av, Nu_av;
#  endif

   nth = g->nth; th = g->th;

/* ============================= Printing result at internal points ========================================== */
   sprintf(resfile, "%02d.%s.%09d.[%E].dat", RUN_CNT, RES_FILE, t_cnt, g->t);

   if((ptr = (FILE *)fopen(resfile, "w")) == (FILE *)0)
     {
      _errno = errno;
      fprintf(g->lptr, "ERROR:: In 'print_result()': Unable to open the result-file '%s' for writing.\n\tSystem errno = %d.\n",\
                  resfile, _errno);
      return -1;
     }

// Printing header for Tecplot.
   fprintf(ptr, "TITLE = \"Re0 = %E, t = %E Sec", g->Re0, g->t);
#  ifdef SOLVE_ENERGY
   if(g->t < t_HEATER_ON) fprintf(ptr, ", (Heater OFF. It will be ON after %E seconds.)", t_HEATER_ON-g->t);
   else fprintf(ptr, ", (Heater ON since %E seconds)", g->t-t_HEATER_ON);
#  endif

   fprintf(ptr, "\"\nVARIABLES = \"r [mm]\", \"z [mm]\", \"th [degree]\", \"u [mm/s]\", \"v [mm/s]\", \"w [mm/s]\", \"Total Velocity [mm/s]\", \"Vorticity_th [rad/sec]\", \"Divergence [kg/m3.sec]\", \"P [Pa]\"");
#  ifdef SOLVE_ENERGY
   fprintf(ptr, ", \"Temperature[ oC]\"");
#  endif
   fprintf(ptr, "\n");

   for(isec = 0; isec < Nsec; isec++)
       {
        u = g->u[isec]; v = g->v[isec]; w = g->w[isec]; p = g->p[isec];
#       ifdef SOLVE_ENERGY
        T = g->T[isec];
#       endif
        r = g->r[isec]; z = g->z[isec]; nr = g->nr[isec]; nz = g->nz[isec];

        fprintf(ptr, "ZONE T=\"Section - %d\", I=%d, J=%d, K=%d, DATAPACKING=POINT\n", isec+1, nr, nz, nth);

        for(j = NthG; j < NthG+nth; j++) for(k = NzG; k < NzG+nz; k++) for(i = NrG; i < NrG+nr; i++)
           {
         /* Evaluating divergence of the velocity field. */
            div = g->density0*DIV_ijk; /* In g/m3.sec */

         /* Evaluating total-velocity. */
            _u = 0.5*(u[i][j][k]+u[i-1][j][k]);
            _v = 0.5*(v[i][j][k]+v[i][j-1][k]);
            _w = 0.5*(w[i][j][k]+w[i][j][k-1]);
            v_total = pow(_u*_u + _v*_v + _w*_w, 0.5);

         /* Evaluating the vorticity in theta-direction. */
            vort_th = VORT_TH_ijk;

            fprintf(ptr, "%E\t%E\t%E\t%E\t%E\t%E\t%E\t%E\t%E\t%E", 0.5*(r[i]+r[i+1])*1000, 0.5*(z[k]+z[k+1])*1000,\
                  0.5*(th[j]+th[j+1])*90/acos(0), _u*1000, _v*1000, _w*1000, v_total*1000, vort_th, div, p[i][j][k]);

#           ifdef SOLVE_ENERGY
            fprintf(ptr, "\t%E", T[i][j][k]-ABS_ZERO);
#           endif

            fprintf(ptr, "\n");
           }
       }
   fclose(ptr);
/* ----------------------------------------------------------------------------------------------------------------------------- */


   if(g->tptr != (FILE *)0) fprintf(g->tptr, "\t'%s' is written!\n", resfile);


/* ======================================= Printing values at the top impingement wall ========================================= */

   sprintf(resfile, "%02d.top.%s.%09d.[%E].dat", RUN_CNT, RES_FILE, t_cnt, g->t);

   if((ptr = (FILE *)fopen(resfile, "w")) == (FILE *)0)
     {
      _errno = errno;
      fprintf(g->lptr, "ERROR:: In 'print_result()': Unable to open the result-file '%s' for writing.\n\tSystem errno = %d.\n",\
                 resfile, _errno);
      return -1;
     }

// Printing header for Tecplot.
   fprintf(ptr, "TITLE = \"Re0 = %E, t = %E Sec", g->Re0, g->t);
#  ifdef SOLVE_ENERGY
   if(g->t < t_HEATER_ON) fprintf(ptr, ", (Heater OFF. It will be ON after %E seconds.)", t_HEATER_ON-g->t);
   else fprintf(ptr, ", (Heater ON since %E seconds)", g->t-t_HEATER_ON);
#  endif

   fprintf(ptr, "\"\nVARIABLES = \"r [mm]\", \"theta [degree]\", \"Radial Stress [N/sq.m]\", \"Tangential Stress [N/sq.m]\"");
#  ifdef SOLVE_ENERGY
   fprintf(ptr, ", \"T(top surface) [oC]\", \"Nu\"");
   T_wall_av = q_av = 0;
#  endif
   fprintf(ptr, "\n");

   u = g->u[1]; v = g->v[1];
#  ifdef SOLVE_ENERGY
   T = g->T[1];
#  endif
   r = g->r[1]; z = g->z[1]; nr = g->nr[1]; nz = g->nz[1];

   fprintf(ptr, "ZONE T=\"Section - %d\", I=%d, J=%d, DATAPACKING=POINT\n", isec+1, nr, nth);

   for(j = NthG; j < NthG + nth; j++) for(i = NrG; i < NrG + nr; i++)
       {
#       ifdef SOLVE_ENERGY
        T_wall = (T[i][j][NzG+nz]+T[i][j][NzG+nz-1])/2;
#       else
        T_wall = T_INLET;
#       endif

       stress_r = kvisc(T_wall)*dens(T_wall)*(u[i][j][NzG+nz-1]+u[i-1][j][NzG+nz-1])/(z[NzG+nz]+z[NzG+nz-1]);
       stress_th = kvisc(T_wall)*dens(T_wall)*(v[i][j][NzG+nz-1]+v[i][j-1][NzG+nz-1])/(z[NzG+nz]+z[NzG+nz-1]);

       fprintf(ptr, "%E\t%E\t%E\t%E", 0.5*(r[i]+r[i+1])*1000, 0.5*(th[j]+th[j+1])*90/acos(0), stress_r, stress_th);

#      ifdef SOLVE_ENERGY
       h = (g->t < t_HEATER_ON ? 0:(i < NrG+g->nr_heater ? TOP_HEATER_HF/(T_wall - T_INLET):0));

       fprintf(ptr, "\t%E\t%E", T_wall - ABS_ZERO, h*2*r[NrG+g->nr_heater]/thm_cond(T_wall));

       T_wall_av += T_wall*0.5*(r[i+1]*r[i+1]-r[i]*r[i])*(th[j+1]-th[j]);
       q_av += (g->t < t_HEATER_ON ? 0:(i < NrG+g->nr_heater ? TOP_HEATER_HF:0))*0.5*(r[i+1]*r[i+1]-r[i]*r[i])*(th[j+1]-th[j]);
#      endif

       fprintf(ptr, "\n");
      }
   fclose(ptr);
/* ----------------------------------------------------------------------------------------------------------------------------- */

   if(g->tptr != (FILE *)0) fprintf(g->tptr, "\t'%s' is written!\n", resfile);


#  ifdef SOLVE_ENERGY
/* Evaluating average value of "Nu". */
   A_total = 0.5*(r[NrG+nr]*r[NrG+nr]-r[NrG]*r[NrG])*(th[NthG+nth]-th[NthG]);
   T_wall_av /= A_total;
   q_av /= A_total;
   h_av = q_av/(T_wall_av - T_INLET);
   Nu_av = h_av*2*r[NrG+g->nr_heater]/thm_cond(T_wall_av);

   if(g->tptr != (FILE *)0) fprintf(g->tptr, "\tNu(av) = %E\n", Nu_av);
#  endif

   return 1;
  }
